Stability of flow through a slowly diverging pipe 



Kirti Chandra Sahu, Rama Govindarajan 
Engineering Mechanics Unit, Jawaharlal Nehru Centre 
for Advanced Scientific Research, Bangalore 560 064, India. 
E-mail: rama@jncasr.ac.in 



February 2, 2008 



Abstract 

Although the critical Reynolds number for linear instability of the laminar flow in a straight pipe 
is infinite, we show that it is finite for a divergent pipe, and approaches infinity as the inverse of the 
divergence angle. The velocity profile at the threshold of inviscid stability is obtained. A non-parallel 
analysis yields linear instability at surprisingly low Reynolds numbers, of about f 50 for a divergence 
of 3 degrees, which would suggest a role for such instabilities in the transition to turbulence. A 
multigrid Poisson equation solver is employed for the basic flow, and an extended eigenvalue method 
for the partial differential equations describing the stability. 



1 Introduction 

The laminar flow through a straight pipe is linearly stable for any Reynolds number [Davev fc Drazinl 
l|l969|) . iLessen et all l)l968(l ]. but, as first demonstrated bv iRevnolds! l)l883(l . transition to turbulence 
usually occurs at a Reynolds number Re, based on the pipe diameter and mean velocity, of around 
2000. By reduc i ng th e external disturbances, it is possible to achieve laminar flow up to Re ~ 10 5 
[Huerre fc Rossil yj)98)] when the pipe is smooth and the flow at the inlet very quiet. The Reynolds 
number u p to which it is p ossible to keep the flow laminar varies inversely as the level of external dis- 
turbance |Hof et all fl20 03)] . Although questions remain about the complete route to turbulence in a 
straight pipe, it s eems likely that the spectrum of linear (stable) modes has a role to play via transient 
algebraic growth [Meseguer fc Trefethenl feonfl : ISchmid fc Hennirigsonl ft 99 A ] of disturbances. It has 
recently been demonstrated both theoretically and experimentally that a nonlinear self-sustaining mech- 
anism leads to the e xistence of travelling waves (and time-periodic sta t es) that appear to pl ay a key role 
in shear turbulence jFaisst fc Eckhardtl 120031 \2004 : iHof eTaH &2QQ4 : IWaleffd (119981120011) 1. 



Our purpose in this paper is to examine the possible role of small local divergences in the transition pro- 
cess. These could have a large effect since linear stability is described by a singular perturbation problem. 
Wheara s a lar ge amount of work has been done on the flow in significantly converging/diverging [e.g. 
iFlorvanl l)2003h ] and in suddenly expanding geometries, we know of no work on the stability of slowly 
diverging pipe flows. Su dden expans i ons ha ve at tracted attentio n beca use of the recirculation zone they 
generate. In particular, iFearn et all l|l990h and ICherdron et all dl978l) study the length o f the recircu- 
lation zone as a function of the Reynolds number, and ISreenivasan fc Strvkowska l|l983j) examine the 
oscillations of the recircul ating bubble and their eff e ct on the flow. Our focus is different, as will become 
clear below. lEaglesl ljl96!Tl) and lEagles fc Weissmanl ft975|) . analysed the Jeffery-Hamel flow generated by 
a slowly diverging channel, and showed by linear parallel stability analysis (the Orr-Sommerfeld equa- 
tion), that the critical Reynolds number decreases by a large amount even for a small divergence angle. 
The divergent pipe is more interesting for several reasons: the critical Reynolds number is infinite for 
an angle of divergence of zero, the Reynolds number is a decreasing function of the streamwise (axial) 
coordinate x, and the flow non-parallelism is larger for a given divergence. In the present work, we 
employ a two-pronged approach. For the mean flow, we derive an axisymmetric Jeffery-Hamel equation 
(AJH), which is valid at small divergence angles. At larger angles of divergence (1° or greater) we solve 
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the Navier-Stokes equations directly in the axisymmetric geometry shown in figure ^ with a divergent 
portion of finite extent. At small angles of divergence and high Reynolds numbers (above 1000), a paral- 
lel flow stability analysis is conducted on the AJH profile, while at lower Reynolds numbers, the partial 
differential equations for non-parallel stability are solved as an extended eigenvalue problem. 

Our main results may be summarised as follows. At low levels of divergence, linear stability is determined 
by a parameter S(x), defined as the product of Re and the slope a of the wall. The basic AJH flow is 
completely described by this parameter. The flow is unstable to the swirl mode for S > 10, so the critical 
Reynolds number approaches infinity as 1/a. At divergences greater than 1°, non-parallel effects are 
found to be quite large, and a non-parallel analysis shows that the flow in a geometry containing a 3° 
divergence is linearly unstable to the swirl mode at Reynolds numbers as low as 150. The following two 
sections describe the basic flow computations and the stability analysis respectively. 



2 The basic flow 

2.1 Axisymmetric Jeffery-Hamel equation 

We begin by noting that unlike in a divergent two-dimensional channel, no similarity flow is possible in a 
divergent pipe. At very low angles of divergence, however, it is possible to derive a one-parameter family 
of velocity profiles, where the parameter 

S = aRe (1) 



varies slowly in the axial coordinate x. Here the slope 

dR{x d ) 

a = ~H 

dxj 



« 1, (2) 



R is the radius of the pipe, and the subscript d stands for a dimensional quantity. Upon eliminating the 
pressure from the axisymmetric momentum equations, we obtain 



v d du d , TT d 2 u d , d 2 u, 



U d ^ " + v d - 



r d dr d dx d dr d dr 2 



d 3 U d 1 dU d , 1 d 2 U d 



Q r 3 r 2 Q rd rd Q r 2 



0, (3) 



where v is the kinematic viscosity, and r is the radial coordinate. The axial and radial velocities, U d and 
V d respectively, may be written in terms of a generalised streamfunction ty d as 

Ud — ^ — , Vd = 5 — , (4) 

rd or d r d ax d 

to satisfy continuity. Substituting this in equation Q, we get 

Zr d* d d 2 v d _ 3 mud%i + r2 ggj d 3 ^ d _ d* d d 2 v d _ r2 d% L d^% L _ 

dxd dr 2 dr d dx d d dr d dr 2 d dx d dr d dr d dxd d dx d dr d 



Td dr\ lTd dr\ + 6rd dr 2 6 dr d 



0. (5) 



The above equation is non-dimensionalised using the local radius R{x d ) and the centreline velocity U c (x d ) 
as scales, e.g., ^ d = U C R 2 ^ . In particular, 

dx d = Rdx. (6) 

The Reynolds number is assumed to be high and the divergence small enough, so that all terms of 
0(Re~ 2 ) or 0(a 2 ) and higher are negligible. The resulting equation in non-dimensional form becomes 

[qV - ra*'] [-3*' + 3r*" - r 2 *'"] - r*' - r*"a - *'a] + r 2 *' 

Uq - 2aW' - ra*'"l - — \r 3 ^ iv - 2r 2 ^"' + 3r^" - 3*'l = 0, (7) 
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Figure 1: Schematic diagram of the divergent pipe used in the numerical simulations, not to scale. 



where 

the Reynolds number is defined as Re(xd) = U c (xd)R(x d )/v, and the primes refer to differentiation with 
respect to r. For the case of near-similar flow, given constant mass flow rate, we may set q = 0, and the 
above equation may be integrated once with respect to r to give 

r 2 U"' = -rU" + (l-4Sr 2 U)U', (8) 

where U = W/r. The boundary conditions are U = at r = 1, and U — 1, U' — at r = 0. Profiles 
obtained from equation (JSJ are compared to those obtained from a numerical axisymmetric Navier-Stokes 
solution in the following subsection. 



2.2 Numerical solution 

The geometry studied here is as shown in the schematic in figure ^ with a straight pipe at the entry, 
followed by an axisymmetric divergent portion, which in turn is followed by a long straight exit portion. 
The length and velocity scales, redefined for convenience in this subsection alone, respectively are the 
radius Ri and the centerline velocity Ui at the inlet. In the case presented here, the divergent portion 
starts at x = 9.4 and ends at x — 91 with a 3° angle of divergence. The total length of the domain is 
L = 120. The axisymmetric Navier-Stokes equations for steady, incompressible, Newtonian flow in the 
streamfunction vorticity formulation, in non-dimensional form, are given by 

« + (t >.V )n .^, (9) 

fi = -V 2 *, (10) 

where Rei = UiRi/v, Q(x,r) is the azimuthal vorticity, U is the velocity vector, and t is time. The 
boundary conditions at the centerline are \& = £1 = V = dU/dr — 0. No-slip and impermeable boundary 
conditions are imposed at the wall. The functional forms of the streamfunction at the centerline, and 
the vorticity at the wall, are described by employing fictitious points outside the domain. At the inlet, a 
parabolic velocity profile is prescribed, while at the outlet the Neumann boundary conditions: d*5f jdx = 0, 
and dQ/dx = are imposed. The reason for using a long exit section, and the consequent increase in 
computational effort, is due to the requirement that the Neumann condition be valid at the exit. 

We begin with a guess solution, where the velocity profile is parabolic at every axial location, and march 
in psuedo-time until a steady-state solution is obtained. The vorticity distribution at each new time 
step is calculated from jnj, adopting a first-order accurate forward differencing in time and second-order 
accurate central differencing in space, on a 34 x 514 grid. The vorticity distribution thus computed is used 
to solve the Poisson equation (|10|l by a Jacobi iterative scheme to obtai n the streamfunc tion everywhere. 
Numerical acceleration is achieved by a six level full-multigrid technique [Fletcherl l|199l[) ] . The procedure 
is repeated until the cumulative change in vorticity during the time step reduces to below = 10~ 10 . 

The axial and radial velocity profiles for Rei = 150 and angle of divergence 3° at different downstream 
locations are compared with those from the AJH profiles in figure EJ^a) and^Jb) respectively. It is seen 
that the AJH profile underpredicts the effect of divergence at x = 22.9 but overpredicts it at x = 46.4. It 
is relevant to mention that the AJH profiles do not correspond to flow through any particular geometry. 
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Figure 2: Comparison of numerically obtained velocity profiles (solid lines) at different axial locations 
(circles: x = 22.9, S = 7.49 and squares: x = 46.4, S — 3.75) with the AJH profiles (dashed lines), (a) 
Axial velocity, (b) Radial velocity. 



At very low angles of divergence, axial variations are slow, the flow attains near-similarity within a short 
downstream distance, and the AJH profiles are expected to predict the real flow well. We find this to be 
the case at angles of divergence less than a degree. The profiles computed here are used in the stability 
calculations, as described in the next section. Incidentally, at higher divergence, regions of separation 
are obtained to very good accuracy by the present method, but are not the subject of discussion here. 
To the contrary, our interest is in finding the smallest divergence at which flow behaviour is completely 
different from that in a straight pipe. 



3 Non-parallel stability analysis 



We now revert to the use of the local radius R(x) and the local centerline velocity U c (x) at a given x as 
scales. Each flow quantity is expressed as the sum of a steady mean and a time-dependent perturbation, 
such as 

u = U(x,r) +u(x,r, 9,t). (11) 

Since the flow under consideration varies significantly in the axial direction, a normal mode form may 
be used only in time and in the azimuthal coordinate 9. In the axial coordinate, the perturbation may 
be expressed as a rapi d ly varying wave-like part scaled by a relatively slowly varying function [see e.g. 
IBertolotti et all l)l992f ): iGovindaraian fc Narasimhal 1 1995h ]. such as 



u, v, w,p] = Real < \u{x, r), v(x, r), w(x, r),p(x, r)] exp 



a(x)dx + n9 — /3dtd 



(12) 



where u, v and w are the axial, radial and the azimuthal velocity perturbations respectively, p is the 
pressure perturbation, a(x) is a local axial wavenumber, n is the number of waves in the azimuthal 
direction, and f3 is the disturbance frequency. Flow quantities in the form (| 1 1|) are substituted in the 
Navier-Stokes and continuity equations, the mean flow equation is subtracted, and nonlinear terms in 
the perturbations are neglected. Since axial variations are slow and the Reynolds number is large, a 
consistent approximation is to retain all terms up to 0(a) and 0(Re~ 1 ) (in the critical and wall layers 
and elsewhere in the pipe) and neglect higher order effects. The result is a set of partial differential 
equations for the perturbation velocities and pressure, each of first order in x and up to second order in 
r, which amounts to a seventh order system in r. These may be expressed in the form 



% r )+^=W(x,r). 



(13) 



Here 4> = [u,v,w,p], (3 = j3dR/U c and the non-zero elements of the 4x4 matrix operators H, Q and B 
are given by 



h n = U 
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m / 1 1 <9 \ <9 

"43 = -5- T + - -5- , = o-; 911 = 322 = 333 = 342 = U, #14 = 1, 

He \r z r or / ar 



and 



fell = &22 = fe33 = fe42 = 1- 



Here U' c = dU c /dx. 



In equatio n i|13|) . w e confirm that if w e set a, U' c and d(f>/dx to zero, we get the parallel stability equations 
of lGiiH {1973) and lLessen et all l)l968|) . The boundary conditions em erge from requiring that all quantities 
vary continuously with r at the centerlinc |P,a,tchelor Rr. Oilll (|1 962h ]. and obey no-slip at the wall: 

u = v = w= p = 0, atr = 0, for 1, (14) 

u = p = 0, v + vw = 0, at r = 0, for n — 1, (15) 

u = v = w = 0, at r = 1. (16) 

Note that for n = 1, we have only six boundary conditions for a seventh order system. We therefore 
generate an extra boundary condition by differentiating the continuity equation with respect to r, and 
using the fact that u(x, 0) = 0, to get 

dv dw 

2— + m— =0 at r = 0, for n = 1. 17 

or or 



Equation (|13|) may be solved as an eigenvalue problem of larger size [Balachandar & Govindarajan, 
preprint, 2005] as described below. The streamwise derivative in equation 1)13(1 couples neighboring axial 
locations in the flow-field to one another. Consider two streamwise locations 1 and 2 separated by an 
incremental distance, i.e., X2 — X\ + Ax. We may write 







O(Ax) 



dx Ax 

Since the dimensional freqency (3d remains constant, (3\ and (3i are related as follows 



^ 2 h i A l^ cl 



We can therefore write (|13(l as 



Hx-Qx/Ax Gil Ax 

-g 2 /Ax n 2 + g 2 /Ax 



01 








kB 2 



(18) 



(19) 



(20) 



(Higher-order approximations to the streamwise derivative could be considered instead of l(18|) and the 
resulting eigensystem would be of correspondingly large size.) The numerical mean flow is interpolated to 
obtain profiles at neighbouring x-locations, with Ace = 0.05. Profiles obtained from computations using 
512 grid points as well as 1024 grid points have been checked to give eigenvaules correct up to 4 decimal 
places. Halving or doubling t he Ax has even less o f an effect on the eigenvalue. Equation (|2lH) is solved by 
a spectral collocation method [Canuto et al\ lll987|) ]. The eigenvalue (3\ is obtained as a complex quantity. 
The complex streamwise wavenumber is iterated until (3\ assumes the desired real value (— (3dR/U c ) at 
a given X\. The axial variation of the wavenumber, da/dx, and the initial guess for a.\ are obtained 



by solving the homogeneous part of the equation ijr3*|) . In subsequent iterations, da/dx is maintained 
constant, since the correction due to the inhomogeneous terms on this quantity is of higher order. It 
is to be noted that the apportionment in (|12f> between the x-depe ndences of a and the ei genfunction is 
arbitrary. There are many ways of performing this apportionment, Bertolotti et and so long 

as the rapid (wavelike) change is included in a ri there is no difference in the prediction of the growth 
of any physical quantity. We have checked that this is the case for the present flow as well. Choosing 
da/dx from the homogeneous part of the equation is one way of including the rapid change into a(x), 
leaving a relatively slow change in u(x). 



We consider downstream growth of disturbances followed at a constant value of the non-dimensional 
radius r. The growth rate g of the nondimensional disturbance kinetic energy, E = l/2(u 2 + v 2 + w 2 ), 
for example, is given by 

1 dE _ 1 dE 



9 E dx 



-2a; + - 

E dx 



(21) 



where E = l/2(uu* + vv* + ww*), the star denotes a complex conjugate. The growth factor for this 
quantity is thus 

E 



E r 



= exp 



g(x)dx 



(22) 



where the subscript cr stands for the critical location, at which g = 0. We see that a disturbance may 
amplify at one r and decay at another. Secondly, one disturbance quantity could be amplifying while 
others decay. 



4 Results and discussion 



The slowest decaying mo d e in a straight pipe is the swirl mode (of azimuthal wave number n = 1) 
|Burridge fc Drazml 1 1969tklCorcos fc Sellarsl (|l9R9() ]. In the diverging pipe, we find that this mode again 
is always the most unstable, all results presented here are for n = 1. We emphasise that for the 3° 
divergence, a non-parallel stability analysis is necessary: with a parallel flow assumption, there is no 
instability until a Reynolds number of about 1000. 



We fir st compare our eigenspectrum for the flow through a straight pipe with that of Sc hmid fc Henningsonl 
l|200lj) . Every eigenvalue matches up to the 10 th decimal place. Next we study the effect of very small 
divergence on the stability behaviour by conducting a parallel flow stability analysis (neglecting non- 
parallel terms in the stability equation) on the AJH profiles. For these parameters, the results have been 
checked to be the same as those from a non-parallel analysis. Figure Eta) shows the critical Reynolds 
number for linear instability as a function of the angle of divergence. At small (but non-zero) divergence 
angle, we find a finite Reynolds number for linear instability. It can be seen that a power-law relationship 
is obeyed. The best fit of the data gives Re cr — 11.2a -0 98 , which is practically indistinguishable from 
Re cr = 10/ a. The critical Reynolds number thus varies as the inverse of the divergence angle. The 
inverse relationship arises because an inviscid mechanism is operational at very high Reynolds numbers, 
and because the AJH velocity profiles JSJ) are described by the product S of Re and a. We thus show 
that there is a limiting velocity profile corresponding to S — 10, and flows where S is greater than 10 
are linearly unstable. When the angle of divergence is less than a degree, the AJH profiles are a good 
approximation of the flow. 

To understand this behaviour better, we write down the equation for inviscid stability (the axisymmetric 
Rayleigh equation) as the divergence goes to zero, by neglecting terms containing Re" 1 and a in l|13|) . 
setting n — 1, and eliminating w, p and u in turn to get 



(U-e) 



3 + a z r 



2J2 









+ 







c?r 2 - 1 
U" - — ^tU' 



r(l + a 2 r 2 ) 



v = 0. 



(23) 



As r — > oo, the above equation, as it should, reduces to the two-dimensional Rayleigh equation. In 
two-dimensional flow, Rayleigh showed that a necessar y cond ition for instability is that U" goes through 
a zero somewhere in the flow [Schmid &: Henning son ( 20 01|) ] . Several improvements on this criterion 
have been made for two-dimensional flow, but there is no proof, as far as we know, of a corresponding 
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Figure 3: (a) Variation of the critical Reynolds number with the tangent a of the divergence angle, at 
small angles of divergence. Symbols: stability analysis; dashed line: best fit; solid line: Re cr = 10/a. (b) 
The inviscid instability function /, in straight pipe flow (S = 0) and for other values of 5*. 




Figure 4: (a) Eigenfunctions u — [u, v, w] and (b) du/dx for Rei = 150, (3d — 0.31 and n = 1 at 
Xd/Ri = 28.1. 



necessary condition for pipe flow. We may however follow a heuristic approach. The quantity / = 
U" — (a 2 r 2 — l)/r/(l + a 2 r 2 )U' is the axisymmetric analogue in 1)23(1 of U" in two-dimensional flow. 
It may therefore be expected that a change of sign in / within the flow would take the flow closer to 
instability. Figure Bib) shows the variation of I with the radial coordinate, I undergoes a sign change if 
S > 2, consistent with expectation. A value of a = 1.26, corresponding to critical instability has been 
used. 

We now examine the behaviour at higher levels of divergence. For the geometry shown in figure^ results 
from a non-parallel spatial stability analysis [equation l|20[)] performed on the numerically obtained profiles 
are presented for Re = 150 at the inlet. The growth rate, as mentioned before, depends on how far the 
monitoring location is from the centerline, and what the quantity being monitored is. An examination 
of equation 1(21(1 shows that the second term on the right hand side determines the r-dependence, and 
comes from the flow quantity under consideration. Typical plots of the eigenfunctions u, v and w, and 
their axial variations, are shown in figures E{a) and (b) respectively. The amplitude of the disturbance 
kinetic energy of the swirl (n = 1) mode is shown in figure El It is seen that the fastest growth occurs 
at r — 0.25, while the disturbance decays near the wall and near the centerline. The Reynolds number 
is a decreasing function of the axial distance, beyond x ~ 50 we find no disturbance that has a positive 
growth rate. At higher Reynolds numbers, there is scope for turbulent flow in the initial portion of the 
divergent section, and relaminarisation downstream. These aspects are being explored experimentally 
(Sood et at, private communication). The sensitivity to small levels of divergence may have implications 
for small scale flows. 

^From this study, we cannot tell at what Reynolds number transition to turbulence will be triggered. 
While an inlet Re of 150 may be too low, it is significant that linear growth has been demonstrated. 
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Figure 5: Amplification of disturbance kinetic energy for Re = 150, n = 1 for typical disturbance 
frequencies, (a) Average across the pipe; (b) r = 0.75; (c) at r = 0.25; (d) at r = 0.07. The axial 
coordinate here is scaled by the inlet radius. 
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At Re — 300 on the other hand, the disturbance kinetic energy grows by a factor of over 60000 so a 
linear mechanism may become important in the transition process. For a continuously diverging pipe, 
since the Reynolds number is a decreasing function of distance, a relaminarisation may occur somewhere 
downstream, so turbulence could be spatially localised. 

In summary, the critical Reynolds number for linear instability of flow in a diverging pipe is finite, and 
can be surprisingly low, with significant disturbance growth rates. The fact that the swirl mode grows 
exponentially indicates that a different route to turbulence from that in a straight pipe is likely. The 
highest growth occurs neither close to the wall nor to th e centreline but in -between. Whether this would 
give rise to structures different from those obtained bv lHof et al. (200^) needs to be investigated. At 
divergences as low as 1 — 2°, the effect of flow non-parallelism is already large, so a non-parallel analysis 
is essential. As the divergence angle goes to zero, the critical Reynolds number approaches infinity as 
X/a. This instability is generated by an inviscid mechanism. 

Grateful thanks are due to Ajay Sood and Narayanan Menon for useful discussions. We thank the Defence 
R&D Organisation, Government of India for financial support. 
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